В любой работе с данными с ростом объемов информации, особенно в больших наборах данных, становится все сложнее находить выбросы и ошибки, которые могут исказить результаты анализа. Исследователи в первую очередь полагаются на визуальный осмотр данных, что удобно для небольших датасетов, однако плохо масштабируется с увеличением количества колонок и наблюдений.
С увеличением объема данных и их сложности стало очевидно, что нужны автоматизированные инструменты. В этом проекте мы рассматриваем наборы инструментов для автоматического поиска выбросов и ошибок в датасетах. Эти инструменты помогут ускорить анализ и сделать его более точным, позволяя исследователям сосредоточиться на более интересных задачах, сделав рутинное ковыряние в данных более эффективным и менее времязатратным.
В качестве примера данных загрузим датасет из проекта 1 (убедитесь, что загрузили эти данные в папку с этим Rmd файлом).
# main_dir <- dirname(rstudioapi::getSourceEditorContext()$path)
# setwd(main_dir)
day_df <- read_csv("day.csv")
## Rows: 731 Columns: 16
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (15): instant, season, yr, mnth, holiday, weekday, workingday, weathers...
## date (1): dteday
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(day_df)
Рассматриваемые подходы и пакеты:
dataMaid: Помогает автоматизировать
проверку качества данных. Он может находить пропуски, ошибки в типах
данных, выбросы, а также автоматически создавать отчёты по результатам
анализа.skimr: Пакет для быстрого суммарного
обзора данных, включая пропущенные значения, распределение переменных и
типы данных.outliers: содержит тесты для оценки
достоверности наличия выбросов.EnvStats: позволяет определить и
оценить потенциальные выбросы.rstatix: позволяет искать многомерные
выбросы на основе расстояния Махаланобиса и скорректированного
расстояния Махаланобиса.boxplot: встроенная функция,
позволяющая строить графики и визуально оценивать наличие выбросов.dataMaid простой пакет, который в 1 клик создаёт rmd-отчёт по вашему датафрейму и конвертирует его в pdf, docx, html.
Большим плюсом является кастомизируемость пакета, в него можно добавить свои функции для автоматического использования, что подробно описано в документации:
The customizability features of
dataMaidessentially comes down to writing suchsummaryFunctions,visualFunctionsandcheckFunctions.
makeDataReport(day_df)
## The default of 'doScale' is FALSE now for stability;
## set options(mc_doScale_quiet=TRUE) to suppress this (once per session) message
## Data report generation is finished. Please wait while your output file is being rendered.
# генерируем отчет в html
makeDataReport(day_df, output = "html", replace = TRUE)
## Data report generation is finished. Please wait while your output file is being rendered.
# подробно о функции
?makeDataReport()
Этот пакет удобен из-за нетрудоёмкости первичного анализа и наглядности генерируемых отчётов, однако не позволяет глубоко закопаться в предоставляемые данные.
Пакет skimr предназначен для предоставления сводной статистики о переменных в дата фреймах, тибблах (tibble), таблицах данных (data tables) и векторах. В базовой библиотеке R наиболее похожую функцию выполняет summary() для векторов и дата фреймов.
Основной функцией пакета skimr является функция skim(), которая предназначена для работы с (группированными) датафреймами, она пытается преобразовать другие объекты в датафреймы, если это возможно. Подобно функции summary(), метод skim() для датафреймов отображает результаты для каждого столбца.
skim(day_df)
| Name | day_df |
| Number of rows | 731 |
| Number of columns | 16 |
| _______________________ | |
| Column type frequency: | |
| Date | 1 |
| numeric | 15 |
| ________________________ | |
| Group variables | None |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dteday | 0 | 1 | 2011-01-01 | 2012-12-31 | 2012-01-01 | 731 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| instant | 0 | 1.00 | 366.00 | 211.17 | 1.00 | 183.50 | 366.00 | 548.50 | 731.00 | ▇▇▇▇▇ |
| season | 0 | 1.00 | 2.50 | 1.11 | 1.00 | 2.00 | 3.00 | 3.00 | 4.00 | ▇▇▁▇▇ |
| yr | 0 | 1.00 | 0.50 | 0.50 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| mnth | 0 | 1.00 | 6.52 | 3.45 | 1.00 | 4.00 | 7.00 | 10.00 | 12.00 | ▇▅▅▅▇ |
| holiday | 0 | 1.00 | 0.03 | 0.17 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| weekday | 0 | 1.00 | 3.00 | 2.00 | 0.00 | 1.00 | 3.00 | 5.00 | 6.00 | ▇▃▃▃▇ |
| workingday | 0 | 1.00 | 0.68 | 0.47 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▃▁▁▁▇ |
| weathersit | 0 | 1.00 | 1.40 | 0.54 | 1.00 | 1.00 | 1.00 | 2.00 | 3.00 | ▇▁▅▁▁ |
| temp | 2 | 1.00 | 0.50 | 0.18 | 0.06 | 0.34 | 0.50 | 0.66 | 0.86 | ▂▇▇▇▅ |
| atemp | 0 | 1.00 | 0.47 | 0.16 | 0.08 | 0.34 | 0.49 | 0.61 | 0.84 | ▂▇▆▇▂ |
| hum | 5 | 0.99 | 0.63 | 0.14 | 0.00 | 0.52 | 0.63 | 0.73 | 0.97 | ▁▁▆▇▂ |
| windspeed | 1 | 1.00 | 0.19 | 0.08 | 0.02 | 0.13 | 0.18 | 0.23 | 0.51 | ▃▇▅▁▁ |
| casual | 0 | 1.00 | 848.18 | 686.62 | 2.00 | 315.50 | 713.00 | 1096.00 | 3410.00 | ▇▆▂▁▁ |
| registered | 1 | 1.00 | 3658.72 | 1559.81 | 20.00 | 2502.25 | 3664.50 | 4783.25 | 6946.00 | ▂▅▇▅▃ |
| cnt | 0 | 1.00 | 4504.35 | 1937.21 | 22.00 | 3152.00 | 4548.00 | 5956.00 | 8714.00 | ▂▅▇▅▃ |
Формат результатов, возвращаемых функцией skim(), представляет собой единый широкий датафрейм, который объединяет результаты анализа и включает в себя несколько атрибутов, а также два важных столбца метаданных:
skim_variable: содержит имена переменных из исходного набора данных. В нашем случае skim_variable — это колонка с именами переменных (Sepal.Length, Sepal.Width и т.д.).
skim_type: обозначает тип данных этих переменных (например, numeric, factor и т.д.). В нашем случае это Variable type: numeric и factor.
Эти столбцы являются ключевыми для объектов класса skim_df и играют важную роль в организации результатов. Если удалить эти столбцы, то объект потеряет свои уникальные свойства и будет преобразован в обычный объект типа tibble. Для проверки, пренадлежит ли объект к классу skim_df, используется функция is_skim_df().
В отличие от функции summary.data.frame(), которая хранит статистику в виде таблицы, объект skim_df имеет свои особенности. Это различие важно, поскольку skim_df можно использовать в пайпах, что облегчает дальнейшую работу с данными. Например, Вы можете выбрать все средние значения переменных или получить сводную статистику для конкретной переменной. Например, получим сводную статистику для столбца Petal.Length:
skim(day_df) %>%
dplyr::filter(skim_variable == "hum") # поддерживается большинство функций из пакета dplyr
| Name | day_df |
| Number of rows | 731 |
| Number of columns | 16 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| hum | 5 | 0.99 | 0.63 | 0.14 | 0 | 0.52 | 0.63 | 0.73 | 0.97 | ▁▁▆▇▂ |
Помимо датафреймов во второй версии пакета skimr функция skim() может анализировать векторы и матрицы путем преобразования их в датафрейм подобно тому как это делает функция as.data.frame():
# Обработка векторов
hum <- day_df$hum # выделяем вектор из данных.
skim(hum)
| Name | hum |
| Number of rows | 731 |
| Number of columns | 1 |
| _______________________ | |
| Column type frequency: | |
| numeric | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| data | 5 | 0.99 | 0.63 | 0.14 | 0 | 0.52 | 0.63 | 0.73 | 0.97 | ▁▁▆▇▂ |
all.equal(skim(hum), skim(as.data.frame(hum)))
## [1] "Attributes: < Component \"df_name\": 1 string mismatch >"
## [2] "Component \"skim_variable\": 1 string mismatch"
# Вывод означает, что объекты, которые мы сравнивали с помощью функции all.equal(), очень похожи, но отличаются в одном атрибуте, конкретно в компоненте df_name.
# Обработка матриц
day_matrix <- as.matrix(day_df)
skim(day_matrix) # аналогична функциям summary.matrix() и colMeans()
| Name | day_matrix |
| Number of rows | 731 |
| Number of columns | 16 |
| _______________________ | |
| Column type frequency: | |
| character | 16 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| instant | 0 | 1.00 | 3 | 3 | 0 | 731 | 0 |
| dteday | 0 | 1.00 | 10 | 10 | 0 | 731 | 0 |
| season | 0 | 1.00 | 1 | 1 | 0 | 4 | 0 |
| yr | 0 | 1.00 | 1 | 1 | 0 | 2 | 0 |
| mnth | 0 | 1.00 | 2 | 2 | 0 | 12 | 0 |
| holiday | 0 | 1.00 | 1 | 1 | 0 | 2 | 0 |
| weekday | 0 | 1.00 | 1 | 1 | 0 | 7 | 0 |
| workingday | 0 | 1.00 | 1 | 1 | 0 | 2 | 0 |
| weathersit | 0 | 1.00 | 1 | 1 | 0 | 3 | 0 |
| temp | 2 | 1.00 | 9 | 9 | 0 | 498 | 0 |
| atemp | 0 | 1.00 | 9 | 9 | 0 | 690 | 0 |
| hum | 5 | 0.99 | 8 | 8 | 0 | 591 | 0 |
| windspeed | 1 | 1.00 | 9 | 9 | 0 | 649 | 0 |
| casual | 0 | 1.00 | 4 | 4 | 0 | 606 | 0 |
| registered | 1 | 1.00 | 4 | 4 | 0 | 678 | 0 |
| cnt | 0 | 1.00 | 4 | 4 | 0 | 696 | 0 |
Пакет skimr так же предоставляет дополнительные функции для работы с данными.
Функция skim_tee() в пакете skimr позволяет получать ту же сводную информацию, что и skim(), но при этом возвращает исходный, неизменённый датафрейм. Это позволяет продолжить работу с оригинальными данными после получения сводной статистики. Посмотрим статистику для day_df, а затем отфильтруем датафрейм по рабочим дням:
day_working <- day_df %>%
skim_tee() %>%
dplyr::filter(workingday == "1")
## ── Data Summary ────────────────────────
## Values
## Name data
## Number of rows 731
## Number of columns 16
## _______________________
## Column type frequency:
## Date 1
## numeric 15
## ________________________
## Group variables None
##
## ── Variable type: Date ─────────────────────────────────────────────────────────
## skim_variable n_missing complete_rate min max median
## 1 dteday 0 1 2011-01-01 2012-12-31 2012-01-01
## n_unique
## 1 731
##
## ── Variable type: numeric ──────────────────────────────────────────────────────
## skim_variable n_missing complete_rate mean sd p0 p25
## 1 instant 0 1 366 211. 1 184.
## 2 season 0 1 2.50 1.11 1 2
## 3 yr 0 1 0.501 0.500 0 0
## 4 mnth 0 1 6.52 3.45 1 4
## 5 holiday 0 1 0.0287 0.167 0 0
## 6 weekday 0 1 3.00 2.00 0 1
## 7 workingday 0 1 0.684 0.465 0 0
## 8 weathersit 0 1 1.40 0.545 1 1
## 9 temp 2 0.997 0.496 0.183 0.0591 0.338
## 10 atemp 0 1 0.474 0.163 0.0791 0.338
## 11 hum 5 0.993 0.628 0.143 0 0.519
## 12 windspeed 1 0.999 0.191 0.0775 0.0224 0.135
## 13 casual 0 1 848. 687. 2 316.
## 14 registered 1 0.999 3659. 1560. 20 2502.
## 15 cnt 0 1 4504. 1937. 22 3152
## p50 p75 p100 hist
## 1 366 548. 731 ▇▇▇▇▇
## 2 3 3 4 ▇▇▁▇▇
## 3 1 1 1 ▇▁▁▁▇
## 4 7 10 12 ▇▅▅▅▇
## 5 0 0 1 ▇▁▁▁▁
## 6 3 5 6 ▇▃▃▃▇
## 7 1 1 1 ▃▁▁▁▇
## 8 1 2 3 ▇▁▅▁▁
## 9 0.5 0.656 0.862 ▂▇▇▇▅
## 10 0.487 0.609 0.841 ▂▇▆▇▂
## 11 0.628 0.730 0.972 ▁▁▆▇▂
## 12 0.181 0.233 0.507 ▃▇▅▁▁
## 13 713 1096 3410 ▇▆▂▁▁
## 14 3664. 4783. 6946 ▂▅▇▅▃
## 15 4548 5956 8714 ▂▅▇▅▃
Функция partition() разделяет результаты, полученные от skim(), на именованный список, где каждый элемент списка соответствует сводной информации для конкретного типа данных:
day_df %>%
skim() %>%
partition()
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dteday | 0 | 1 | 2011-01-01 | 2012-12-31 | 2012-01-01 | 731 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| instant | 0 | 1.00 | 366.00 | 211.17 | 1.00 | 183.50 | 366.00 | 548.50 | 731.00 | ▇▇▇▇▇ |
| season | 0 | 1.00 | 2.50 | 1.11 | 1.00 | 2.00 | 3.00 | 3.00 | 4.00 | ▇▇▁▇▇ |
| yr | 0 | 1.00 | 0.50 | 0.50 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| mnth | 0 | 1.00 | 6.52 | 3.45 | 1.00 | 4.00 | 7.00 | 10.00 | 12.00 | ▇▅▅▅▇ |
| holiday | 0 | 1.00 | 0.03 | 0.17 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| weekday | 0 | 1.00 | 3.00 | 2.00 | 0.00 | 1.00 | 3.00 | 5.00 | 6.00 | ▇▃▃▃▇ |
| workingday | 0 | 1.00 | 0.68 | 0.47 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▃▁▁▁▇ |
| weathersit | 0 | 1.00 | 1.40 | 0.54 | 1.00 | 1.00 | 1.00 | 2.00 | 3.00 | ▇▁▅▁▁ |
| temp | 2 | 1.00 | 0.50 | 0.18 | 0.06 | 0.34 | 0.50 | 0.66 | 0.86 | ▂▇▇▇▅ |
| atemp | 0 | 1.00 | 0.47 | 0.16 | 0.08 | 0.34 | 0.49 | 0.61 | 0.84 | ▂▇▆▇▂ |
| hum | 5 | 0.99 | 0.63 | 0.14 | 0.00 | 0.52 | 0.63 | 0.73 | 0.97 | ▁▁▆▇▂ |
| windspeed | 1 | 1.00 | 0.19 | 0.08 | 0.02 | 0.13 | 0.18 | 0.23 | 0.51 | ▃▇▅▁▁ |
| casual | 0 | 1.00 | 848.18 | 686.62 | 2.00 | 315.50 | 713.00 | 1096.00 | 3410.00 | ▇▆▂▁▁ |
| registered | 1 | 1.00 | 3658.72 | 1559.81 | 20.00 | 2502.25 | 3664.50 | 4783.25 | 6946.00 | ▂▅▇▅▃ |
| cnt | 0 | 1.00 | 4504.35 | 1937.21 | 22.00 | 3152.00 | 4548.00 | 5956.00 | 8714.00 | ▂▅▇▅▃ |
Функция yank() в свою очередь извлекает подтаблицу, соответствующую конкретному типу данных. Воспринимайте это как о dplyr::select для типов столбцов в исходных данных:
day_df %>%
skim() %>%
yank("numeric")
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| instant | 0 | 1.00 | 366.00 | 211.17 | 1.00 | 183.50 | 366.00 | 548.50 | 731.00 | ▇▇▇▇▇ |
| season | 0 | 1.00 | 2.50 | 1.11 | 1.00 | 2.00 | 3.00 | 3.00 | 4.00 | ▇▇▁▇▇ |
| yr | 0 | 1.00 | 0.50 | 0.50 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| mnth | 0 | 1.00 | 6.52 | 3.45 | 1.00 | 4.00 | 7.00 | 10.00 | 12.00 | ▇▅▅▅▇ |
| holiday | 0 | 1.00 | 0.03 | 0.17 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| weekday | 0 | 1.00 | 3.00 | 2.00 | 0.00 | 1.00 | 3.00 | 5.00 | 6.00 | ▇▃▃▃▇ |
| workingday | 0 | 1.00 | 0.68 | 0.47 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▃▁▁▁▇ |
| weathersit | 0 | 1.00 | 1.40 | 0.54 | 1.00 | 1.00 | 1.00 | 2.00 | 3.00 | ▇▁▅▁▁ |
| temp | 2 | 1.00 | 0.50 | 0.18 | 0.06 | 0.34 | 0.50 | 0.66 | 0.86 | ▂▇▇▇▅ |
| atemp | 0 | 1.00 | 0.47 | 0.16 | 0.08 | 0.34 | 0.49 | 0.61 | 0.84 | ▂▇▆▇▂ |
| hum | 5 | 0.99 | 0.63 | 0.14 | 0.00 | 0.52 | 0.63 | 0.73 | 0.97 | ▁▁▆▇▂ |
| windspeed | 1 | 1.00 | 0.19 | 0.08 | 0.02 | 0.13 | 0.18 | 0.23 | 0.51 | ▃▇▅▁▁ |
| casual | 0 | 1.00 | 848.18 | 686.62 | 2.00 | 315.50 | 713.00 | 1096.00 | 3410.00 | ▇▆▂▁▁ |
| registered | 1 | 1.00 | 3658.72 | 1559.81 | 20.00 | 2502.25 | 3664.50 | 4783.25 | 6946.00 | ▂▅▇▅▃ |
| cnt | 0 | 1.00 | 4504.35 | 1937.21 | 22.00 | 3152.00 | 4548.00 | 5956.00 | 8714.00 | ▂▅▇▅▃ |
В документации представлены и другие функции из этого пакета, выше перечислены основные.
Отображение результатов skim()
Объект skim_df, получаемый от skim(), представлен в виде широкого датафрейма. По умолчанию его отображение осуществляется с помощью функции print.skim_df(). Для документов, создаваемых с использованием knitr, пакет skimr предоставляет специальный метод knit_print, который автоматически форматирует вывод. Чтобы использовать этот метод, последний оператор в вашем коде должен возвращать объект skim_df.
skim(day_df)
| Name | day_df |
| Number of rows | 731 |
| Number of columns | 16 |
| _______________________ | |
| Column type frequency: | |
| Date | 1 |
| numeric | 15 |
| ________________________ | |
| Group variables | None |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| dteday | 0 | 1 | 2011-01-01 | 2012-12-31 | 2012-01-01 | 731 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| instant | 0 | 1.00 | 366.00 | 211.17 | 1.00 | 183.50 | 366.00 | 548.50 | 731.00 | ▇▇▇▇▇ |
| season | 0 | 1.00 | 2.50 | 1.11 | 1.00 | 2.00 | 3.00 | 3.00 | 4.00 | ▇▇▁▇▇ |
| yr | 0 | 1.00 | 0.50 | 0.50 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▇▁▁▁▇ |
| mnth | 0 | 1.00 | 6.52 | 3.45 | 1.00 | 4.00 | 7.00 | 10.00 | 12.00 | ▇▅▅▅▇ |
| holiday | 0 | 1.00 | 0.03 | 0.17 | 0.00 | 0.00 | 0.00 | 0.00 | 1.00 | ▇▁▁▁▁ |
| weekday | 0 | 1.00 | 3.00 | 2.00 | 0.00 | 1.00 | 3.00 | 5.00 | 6.00 | ▇▃▃▃▇ |
| workingday | 0 | 1.00 | 0.68 | 0.47 | 0.00 | 0.00 | 1.00 | 1.00 | 1.00 | ▃▁▁▁▇ |
| weathersit | 0 | 1.00 | 1.40 | 0.54 | 1.00 | 1.00 | 1.00 | 2.00 | 3.00 | ▇▁▅▁▁ |
| temp | 2 | 1.00 | 0.50 | 0.18 | 0.06 | 0.34 | 0.50 | 0.66 | 0.86 | ▂▇▇▇▅ |
| atemp | 0 | 1.00 | 0.47 | 0.16 | 0.08 | 0.34 | 0.49 | 0.61 | 0.84 | ▂▇▆▇▂ |
| hum | 5 | 0.99 | 0.63 | 0.14 | 0.00 | 0.52 | 0.63 | 0.73 | 0.97 | ▁▁▆▇▂ |
| windspeed | 1 | 1.00 | 0.19 | 0.08 | 0.02 | 0.13 | 0.18 | 0.23 | 0.51 | ▃▇▅▁▁ |
| casual | 0 | 1.00 | 848.18 | 686.62 | 2.00 | 315.50 | 713.00 | 1096.00 | 3410.00 | ▇▆▂▁▁ |
| registered | 1 | 1.00 | 3658.72 | 1559.81 | 20.00 | 2502.25 | 3664.50 | 4783.25 | 6946.00 | ▂▅▇▅▃ |
| cnt | 0 | 1.00 | 4504.35 | 1937.21 | 22.00 | 3152.00 | 4548.00 | 5956.00 | 8714.00 | ▂▅▇▅▃ |
Пакет outliers содержит в себе несколько статистических тестов, позволяющих оценить наличие выбросов в принципе или задать их вручную и оценить достоверность этих выбросов.
позволяет оценить, является ли максимальное или минимальное значение выбросом и оценивает уровень значимости. Работает на выборках размером больше \(6\).
По умолчанию считает нижний порог, чтобы посчитать статистику
верхнего, нужно добавить параметр opposite = TRUE.
day_df = read.csv('day.csv')
grubbs.test(day_df$hum)
##
## Grubbs test for one outlier
##
## data: day_df$hum
## G = 4.39737, U = 0.97329, p-value = 0.003491
## alternative hypothesis: lowest value 0 is an outlier
grubbs.test(day_df$hum, opposite = TRUE)
##
## Grubbs test for one outlier
##
## data: day_df$hum
## G = 2.41300, U = 0.99196, p-value = 1
## alternative hypothesis: highest value 0.9725 is an outlier
В данном случае мы отвергаем, что значение \(0.9725\) это выброс, но принимаем, что \(0\) это выброс. Смотрим на графике:
out_day_df <- boxplot.stats(day_df$hum)$out
boxplot(day_df$hum)
mtext(paste("Outlier: ", paste(out_day_df, collapse = ", ")))
в отличие от теста Грабса, работает на маленьких выборках \(<=25\).
subday_df <- day_df[1:25, ] # отбираем из датасета первых 20 значений.
dixon.test(subday_df$hum)
##
## Dixon test for outliers
##
## data: subday_df$hum
## Q = 0.28209, p-value = 0.558
## alternative hypothesis: highest value 0.861667 is an outlier
В данном случае p.value слишком высок, значит мы отвергаем альтернативную гипотезу о том, что значение из нашей подвыборки является выбросом.
Можно перепроверить это с помощью boxplot.
out <- boxplot.stats(subday_df$hum)$out
boxplot(subday_df$hum)
mtext(paste("Outlier: ", paste(out, collapse = ", ")))
Пакет EnvStats содержит тест Рознера, который, в отличие от двух предыдущих тестов, используется для поиска сразу нескольких выбросов. Работает на выборках больше \(20\).
Здесь мы можем протестировать сразу несколько значений. Тест покажет, являются ли выбросы
Посмотрим на boxplot, чтобы определить потенциальное количество выбросов.
out_day_df <- boxplot.stats(day_df$hum)$out
boxplot(day_df$hum)
mtext(paste("Outlier: ", paste(out_day_df, collapse = ", ")))
При применении теста мы видим, что реальным выбросом является только 0. Второе значение не является выбросом. Мы можем узнать номер выброса в таблице при выводе теста, а вывод all.stats показывает нам произведённые расчёты. Функция автоматически удаляет NA из данных перед подсчётом квантилей.
test <- rosnerTest(day_df$hum, k = 2)
## Warning in rosnerTest(day_df$hum, k = 2): 5 observations with NA/NaN/Inf in 'x'
## removed.
test$all.stats
Пакет rstatix содержит функции для выявления одномерных и многомерных выбросов. Ищет с помощью boxplot и квартилей. Определяет выбросы и экстремальные выбросы с помощью квантилей (выше Q3 + 3 x IQR или ниже Q1 - 3 x IQR).
out_rstatix = identify_outliers(data = day_df, variable = 'hum')
out_rstatix
Пакет MVN ищет многомерные выбросы на основе расстояния Махаланобиса и скорректированного расстояния Махаланобиса.
С помощью функции mvn из пакета можно нарисовать графики Q-Q plot с выбросами на них по каждому из наблюдений, а также скорректированный график \(\chi\)-square.
# install.packages("MVN")
library(MVN)
sub_days = day_df[c(7, 10, 11, 12)]
result = mvn(data = sub_days, subset = "weekday", mvnTest = "hz",
univariateTest = "AD", univariatePlot = "histogram",
multivariatePlot = "qq", multivariateOutlierMethod = "adj",
showOutliers = TRUE, showNewData = TRUE)
result$multivariateOutliers
## $`0`
## Observation Mahalanobis Distance Outlier
## 408 408 95.516 TRUE
## 9 9 82.867 TRUE
## 23 23 49.306 TRUE
## 730 730 45.222 TRUE
## 380 380 41.246 TRUE
## 387 387 33.693 TRUE
## 429 429 27.186 TRUE
## 555 555 25.030 TRUE
## 415 415 24.226 TRUE
## 394 394 24.115 TRUE
## 205 205 23.205 TRUE
## 422 422 22.732 TRUE
## 51 51 18.688 TRUE
## 86 86 17.486 TRUE
## 16 16 16.551 TRUE
## 219 219 14.922 TRUE
## 268 268 14.211 TRUE
## 401 401 12.809 TRUE
## 352 352 10.556 TRUE
##
## $`1`
## Observation Mahalanobis Distance Outlier
## 10 10 82.805 TRUE
## 66 66 81.149 TRUE
## 367 367 75.787 TRUE
## 17 17 66.147 TRUE
## 24 24 62.073 TRUE
## 381 381 60.689 TRUE
## 87 87 54.167 TRUE
## 31 31 52.572 TRUE
## 52 52 50.578 TRUE
## 395 395 46.695 TRUE
## 430 430 44.501 TRUE
## 479 479 44.002 TRUE
## 409 409 42.973 TRUE
## 416 416 38.610 TRUE
## 731 731 35.834 TRUE
## 675 675 32.887 TRUE
## 206 206 32.539 TRUE
## 360 360 24.724 TRUE
## 192 192 21.450 TRUE
## 353 353 20.927 TRUE
## 45 45 18.502 TRUE
## 423 423 16.820 TRUE
## 388 388 14.583 TRUE
## 563 563 13.999 TRUE
## 374 374 13.338 TRUE
## 80 80 11.888 TRUE
## 269 269 11.480 TRUE
## 73 73 11.414 TRUE
## 402 402 10.818 TRUE
## 549 549 9.983 TRUE
##
## $`2`
## Observation Mahalanobis Distance Outlier
## 368 368 46.817 TRUE
## 39 39 28.605 TRUE
## 46 46 20.559 TRUE
## 88 88 15.249 TRUE
## 452 452 13.600 TRUE
## 53 53 12.317 TRUE
## 200 200 11.144 TRUE
##
## $`3`
## Observation Mahalanobis Distance Outlier
## 383 383 30.201 TRUE
## 369 369 29.074 TRUE
## 726 726 28.533 TRUE
## 26 26 24.137 TRUE
## 40 40 23.145 TRUE
## 677 677 21.085 TRUE
## 201 201 20.319 TRUE
## 362 362 20.022 TRUE
## 5 5 16.828 TRUE
## 208 208 14.637 TRUE
## 61 61 13.635 TRUE
## 222 222 12.921 TRUE
## 684 684 11.593 TRUE
## 698 698 11.515 TRUE
## 614 614 11.344 TRUE
## 33 33 11.211 TRUE
## 334 334 9.887 TRUE
## 54 54 9.813 TRUE
##
## $`4`
## Observation Mahalanobis Distance Outlier
## 202 202 110.358 TRUE
## 13 13 48.368 TRUE
## 727 727 47.870 TRUE
## 251 251 46.992 TRUE
## 265 265 37.518 TRUE
## 34 34 34.559 TRUE
## 69 69 29.710 TRUE
## 90 90 28.873 TRUE
## 41 41 25.555 TRUE
## 83 83 24.289 TRUE
## 62 62 22.232 TRUE
## 384 384 21.661 TRUE
## 342 342 17.089 TRUE
## 153 153 16.357 TRUE
## 20 20 15.495 TRUE
## 321 321 13.871 TRUE
## 405 405 11.745 TRUE
## 55 55 11.434 TRUE
## 545 545 11.323 TRUE
## 678 678 10.811 TRUE
##
## $`5`
## Observation Mahalanobis Distance Outlier
## 595 595 1914.826 TRUE
## 266 266 45.125 TRUE
## 203 203 40.278 TRUE
## 21 21 21.654 TRUE
## 252 252 19.709 TRUE
## 378 378 13.952 TRUE
##
## $`6`
## Observation Mahalanobis Distance Outlier
## 302 302 32.619 TRUE
## 421 421 32.106 TRUE
## 722 722 29.848 TRUE
## 694 694 27.903 TRUE
## 407 407 20.640 TRUE
## 8 8 19.133 TRUE
## 22 22 18.298 TRUE
## 386 386 17.634 TRUE
## 379 379 13.414 TRUE
## 204 204 12.757 TRUE
## 85 85 11.826 TRUE
## 435 435 11.407 TRUE
## 540 540 11.367 TRUE
## 351 351 11.349 TRUE
## 50 50 9.929 TRUE
Попробуем “полуавтоматический” поиск выбросов. Такой способ позволяет контролировать каждый этап. Воспользуемся встроенным методом boxplot.
С помощью Boxplot можно не только увидеть выбросы графически, но и сохранить их для последующего отсечения из данных. Для визуализации возьмём рандомные данные.
# создаём данные: с выбросами и без.
dirty_data <- c(4.358015, 4.489513, 4.560919, 4.613810, 4.599738, 4.621614, 4.633119, 4.616862, 4.754681, 4.849953, 4.945791, 5.019631, 4.805033, 4.989170, 5.024305, 5.065325, 4.970247, 4.998086, 5.096887, 4.977657, 4.888269, 3.479053, 2.878145)
data <- c(2.633213, 2.654674, 2.746650, 2.657763, 2.525229, 2.549804, 2.537088, 1.974909, 1.838017, 1.791683, 1.782088, 1.664908, 1.689402, 1.688826, 1.661763, 1.734322, 1.744875, 1.710471, 1.735690, 1.800677, 1.607354, 1.896810, 2.294757)
plot(data, dirty_data)
boxplot(data)
boxplot(dirty_data)
Визуально можем найти выбросы на графике. Чтобы удалить их, можно воспользоваться
boxplot не только рисует картинку, но и сохраняет все ее параметры в объекте, из которого мы можем их достать. Выбросы хранятся тут:
boxplot.stats(y)$out
Можно найти индексы точек выбросов.
ind <- which(dirty_data %in% boxplot.stats(dirty_data)$out)
ind
## [1] 22 23
Сохраним координаты точек выбросов в отдельном dataframe
out <- data.frame(data=data[ind], dirty_data=dirty_data[ind])
out
Если создать переменную с графиком (одним или несколькими), при выводе переменной мы получим следующие данные:
res <- boxplot(dirty_data)
res
## $stats
## [,1]
## [1,] 4.358015
## [2,] 4.606774
## [3,] 4.805033
## [4,] 4.983413
## [5,] 5.096887
##
## $n
## [1] 23
##
## $conf
## [,1]
## [1,] 4.680948
## [2,] 4.929118
##
## $out
## [1] 3.479053 2.878145
##
## $group
## [1] 1 1
##
## $names
## [1] "1"
И с помощью созданной переменной можно снова создать boxplot.
bxp(res)
Теперь визуализируем их на обычном графике:
plot(data, dirty_data, col='green', pch=18, ylim=c(0,max(dirty_data)))
points(out$data, out$dirty_data, col='red', pch=18)
Удалим точки из данных.
clean_data <- dirty_data[-ind]
data_data <- data[-ind]
clean_data
## [1] 4.358015 4.489513 4.560919 4.613810 4.599738 4.621614 4.633119 4.616862
## [9] 4.754681 4.849953 4.945791 5.019631 4.805033 4.989170 5.024305 5.065325
## [17] 4.970247 4.998086 5.096887 4.977657 4.888269
data_data
## [1] 2.633213 2.654674 2.746650 2.657763 2.525229 2.549804 2.537088 1.974909
## [9] 1.838017 1.791683 1.782088 1.664908 1.689402 1.688826 1.661763 1.734322
## [17] 1.744875 1.710471 1.735690 1.800677 1.607354
Сравним боксплоты с выбросами и без них:
boxplot(dirty_data, data_data)
boxplot(clean_data, data_data)
Сравним графики с выбросами и без них:
plot(data, dirty_data, col='green', pch=18, ylim=c(0,max(dirty_data)))
points(out$data, out$dirty_data, col='red', pch=18)
plot(data_data, clean_data, col='green', pch=18, ylim = c(0, max(clean_data)))
заключаем что вы биоинформатик (^-人-^)